Antenna reader example

from radiocalibrationtoolkit import *
from healpy.newvisufunc import projview
import plotly.graph_objects as go

# This ensures Plotly output works in multiple places:
# plotly_mimetype: VS Code notebook UI
# notebook: "Jupyter: Export to HTML" command in VS Code
# See https://plotly.com/python/renderers/#multiple-renderers
import plotly.io as pio
pio.renderers.default = "plotly_mimetype+notebook"

[INFO] LFmap: Import successful.

layout_settings = dict(
    xaxis=dict(title="<b>azimuth</b>", tickprefix="<b>", ticksuffix="</b>", dtick=30),
    yaxis=dict(
        title="<b>zenith angle</b>",
        tickprefix="<b>",
        ticksuffix="</b>",
        range=(0, 90),
        tick0=0,
        dtick=10,
        autorange=False,
    ),
    coloraxis=dict(
        colorbar=dict(
            title=dict(
                text="<b>VEL</b>",
                side="right",
            ),
            tickprefix="<b>",
            ticksuffix="</b>",
        ),
    ),
    font=dict(
        size=15,
        color="black",
    ),
)

# create antenna instance
antenna_inst = AntennaPattern("./antenna_setup_files/SALLA_EW.xml")

[INFO] Your keys are: ['EAHTheta_amp', 'EAHTheta_phase', 'EAHPhi_amp', 'EAHPhi_phase', 'absolute'] Use them as the 'quantity'

# XML values in dictionary
# antenna_inst.get_raw()

# get antenna gain for unpolarized emission
df = antenna_inst.get(frequency=45, quantity="absolute")

# plot
fig = px.imshow(df.T.iloc[::-1, :], width=600, aspect="cube")
fig.update_layout(**layout_settings)
fig.show()

# interpolate
df = antenna_inst.get(frequency=45, quantity='absolute', interp_phi=np.linspace(0, 360, 200), interp_theta=np.linspace(0,90, 100))

# plot
fig = px.imshow(df.T.iloc[::-1, :], width=600, aspect="cube")
fig.update_layout(**layout_settings)
fig.show()

## show antenna pattern as volumetric data with slices in all axis

# quantity = "absolute"
# quantity = "EAHPhi_amp"
quantity = "EAHTheta_amp"

def create_volume(antenna_inst, quantity = "EAHTheta_amp", freq_range=np.arange(30, 81, 3)):
    """
    Create a volume based on the given quantity and frequency range.

    Parameters
    ----------
    antenna_inst : AntennaPattern
        The antenna pattern instance to use for retrieving the data.
    quantity : str, optional
        The quantity to use. Valid input strings are "absolute", "EAHPhi_amp", "EAHTheta_amp".
        Default is "EAHTheta_amp".
    freq_range : np.ndarray, optional
        The range of frequencies. Default is np.arange(30, 81, 3).

    Returns
    -------
    tuple
        A tuple containing PHI, FREQ, THETA, and volume arrays.
    """

    df = antenna_inst.get(frequency=30, quantity=quantity)
    if quantity == 'absolute':
        quantity_label = "<|H|>"
    elif "Phi_amp" in quantity:
        quantity_label = "|H<sub>Φ</sub>|"
    elif "Theta_amp" in quantity:
        quantity_label = "|H<sub>ϴ</sub>|"

    values = np.array([])
    volume_data = []
    for f in freq_range:
        df = antenna_inst.get(frequency=f, quantity=quantity)
        values = np.append(values, df.values.flatten())
        volume_data.append(df.values)

    volume_data = np.asarray(volume_data)
    PHI, FREQ, THETA = np.meshgrid(df.index.values, freq_range, df.columns.values)
    return PHI, FREQ, THETA, volume_data, quantity_label


# some settings
font = "Arial Black"
colorscale='jet'

scene = dict(
    xaxis=dict(
        title="<b>Zenith angle [°]</b>",
        color="black",
        dtick="30"
    ),
    yaxis=dict(
        title="<b>Azimuth [°]</b>",
        color="black",
        dtick="90"
    ),
    zaxis=dict(
        title="<b>frequency [MHz]</b>",
        tickvals=[40, 50, 60, 70, 80],
    ),
    aspectratio=dict(x=0.75, y=0.75,
                     z=0.75),  # Adjust the aspect ratio as needed
)

def quantity_label2colorbar_dict(quantity_label):
    return dict(
                len=0.75,
                title=dict(
                    text="<b>" + quantity_label + " [m]</b>",
                    side="right",
                ),
                tickprefix="<b>",
                ticksuffix="</b>",
            )


def create_volume_trace(
        opacity=1,
        surface_fill=0.001,
        surface_count=1,
        opacityscale="uniform",
        caps=dict(x_show=False, y_show=False, z_show=False),
):
    """
    Create a go.Volume trace based on the given parameters.

    Parameters
    ----------
    opacity : float, optional
        The opacity of the volume. Default is 1.
    surface_fill : float, optional
        The surface fill value. Default is 0.001.
    surface_count : int, optional
        The surface count value. Default is 1.
    opacityscale : str, optional
        The opacity scale. Default is "uniform".
    caps : dict, optional
        The caps configuration. Default is dict(x_show=False, y_show=False, z_show=False).

    Returns
    -------
    go.Volume
        A go.Volume trace object.
    """
    return go.Volume(
        y=PHI.flatten(),
        x=THETA.flatten(),
        z=FREQ.flatten(),
        value=np.around(volume_data.flatten(), 3),
        cmin=0,
        cmax=2,
        opacity=opacity,
        # isomax=1.3,
        surface_fill=surface_fill,
        surface_count=surface_count,
        colorscale="jet",
        opacityscale=opacityscale,
        caps=caps,
        colorbar=quantity_label2colorbar_dict(quantity_label),
    )


# create volume data
PHI, FREQ, THETA, volume_data, quantity_label = create_volume(antenna_inst, quantity = "EAHTheta_amp")

# define some slices
slice_types = {
    "no_slices": [],  # Neutral entry for no slices
    "slices_z": [45, 60],
    "slices_x": [30, 65],
    "slices_y": [0, 270],
}

# make plots
for slice_type, locations in slice_types.items():
    if slice_type == "no_slices":
        fig = go.Figure()
        fig.add_trace(
            create_volume_trace(opacity=0.8,
                                surface_fill=1,
                                opacityscale="max",
                                surface_count=40,
                                caps=dict(x_show=True,
                                          y_show=True,
                                          z_show=True)))
    else:
        fig = go.Figure()
        slice_trace = create_volume_trace()
        slice_trace.update(
            **{slice_type: dict(show=True, locations=locations)})
        fig.add_trace(slice_trace)

    fig.update_layout(
        scene=scene,
        margin=dict(r=50, b=10, l=10, t=10),
        height=600,
        autosize=False,
        font=dict(family=font, size=18, color="black"),
    )
    fig.show()

## show antenna pattern as volumetric data with slices in all axis
# version of plots with sliders

def show_slices(x, y, z, volume_data, slice_type='x'):
    """
    Show slices of volumetric data using the given parameters.

    Parameters
    ----------
    x : np.ndarray
        The x-coordinate values.
    y : np.ndarray
        The y-coordinate values.
    z : np.ndarray
        The z-coordinate values.
    volume_data : np.ndarray
        The volumetric data.
    slice_type : str, optional
        The type of slice. Valid values are 'x', 'y', or 'z'. Default is 'x'.

    Returns
    -------
    None
        This function does not return anything. It displays the plot.

    """
    x_delta = np.diff(x)[0]
    y_delta = np.diff(y)[0]
    z_delta = np.diff(z)[0]

    cmin = 0
    cmax = 2
    colorscale = 'jet'

    if slice_type == 'x':
        slider_label = x
        fig = go.Figure(frames=[
            go.Frame(
                data=go.Surface(x=np.ones(x.size) * x[k],
                                y=y,
                                z=(np.ones((z.size, y.size)) *
                                   z[:, np.newaxis]).T,
                                surfacecolor=(volume_data[:, :, k]).T,
                                cmin=cmin,
                                cmax=cmax),
                name=str(
                    k
                )
            ) for k in range(x.size)
        ])

        # default fig
        fig.add_trace(
            go.Surface(
                x=np.ones(x.size) * x[0],
                y=y,
                z=(np.ones((z.size, y.size)) *
                   z[:, np.newaxis]).T,
                surfacecolor=(volume_data[:, :, 0]).T,
                colorscale=colorscale,
                cmin=cmin,
                cmax=cmax,
                colorbar=quantity_label2colorbar_dict(quantity_label),
            ))

    elif slice_type == 'y':
        slider_label = y
        fig = go.Figure(frames=[
            go.Frame(
                data=go.Surface(
                    x=x,
                    y=np.ones(y.size) * y[k],
                    z=np.ones((z.size, x.size)) *
                    z[:, np.newaxis],
                    surfacecolor=(volume_data[:, y.size - 1 - k, :]),
                    cmin=cmin,
                    cmax=cmax),
                name=str(
                    k
                )
            ) for k in range(y.size)
        ])

        # default fig
        fig.add_trace(
            go.Surface(
                x=x,
                y=np.ones(y.size) * y[0],
                z=np.ones((z.size, x.size)) *
                z[:, np.newaxis],
                surfacecolor=(volume_data[:, 0, :]),
                colorscale=colorscale,
                cmin=cmin,
                cmax=cmax,
                colorbar=quantity_label2colorbar_dict(quantity_label),
            ))
    elif slice_type == 'z':
        slider_label = z

        fig = go.Figure(frames=[
            go.Frame(
                data=go.Surface(x=x,
                                y=y,
                                z=z[k] * np.ones((y.size, x.size)),
                                surfacecolor=np.flipud(volume_data[z.size - 1 - k]),
                                cmin=cmin,
                                cmax=cmax),
                name=str(
                    k
                )
            ) for k in range(z.size)
        ])

        # default fig
        fig.add_trace(
            go.Surface(
                x=x,
                y=y,
                z=z[0] * np.ones((y.size, x.size)),
                surfacecolor=np.flipud(volume_data[z.size - 1]),
                colorscale=colorscale,
                cmin=cmin,
                cmax=cmax,
                colorbar=quantity_label2colorbar_dict(quantity_label),
            ))


    def frame_args(duration):
        return {
            "frame": {
                "duration": duration
            },
            "mode": "immediate",
            "fromcurrent": True,
            "transition": {
                "duration": duration,
                "easing": "linear"
            },
        }


    sliders = [{
        "pad": {
            "b": 10,
            "t": 60
        },
        "len":
        0.9,
        "x":
        0.1,
        "y":
        0,
        "steps": [{
            "args": [[f.name], frame_args(0)],
            "label": str(slider_label[k]),
            "method": "animate",
        } for k, f in enumerate(fig.frames)],
    }]

    # Layout
    fig.update_layout(
        title='<br>Slices in volumetric data',
        width=800,
        height=800,
        scene=dict(
            yaxis=dict(range=[y[0] - y_delta, y[-1] + y_delta],
                       autorange=False),  # Set the y-axis range from -5 to +5
            zaxis=dict(range=[z[0] - z_delta, z[-1] + z_delta], autorange=False),
            xaxis=dict(range=[x[0] - x_delta, x[-1] + x_delta], autorange=False),
            aspectratio=dict(x=1, y=1, z=1),
        ),
        updatemenus=[{
            "buttons": [
                {
                    "args": [None, frame_args(50)],
                    "label": "&#9654;",  # play symbol
                    "method": "animate",
                },
                {
                    "args": [[None], frame_args(0)],
                    "label": "&#9724;",  # pause symbol
                    "method": "animate",
                },
            ],
            "direction":
            "left",
            "pad": {
                "r": 10,
                "t": 70
            },
            "type":
            "buttons",
            "x":
            0.1,
            "y":
            0,
        }],
        sliders=sliders)

    fig.update_layout(
        scene=scene,
        margin=dict(r=50, b=10, l=10, t=10),
        font=dict(family=font, size=18, color="black"),
    )

    fig.show()

# create slices with sliders
z = FREQ[:,0,0]
y = PHI[0,:,0]
x = THETA[0,0,:]
PHI, FREQ, THETA, volume_data, quantity_label = create_volume(antenna_inst, quantity = "EAHTheta_amp")
show_slices(x, y, z, volume_data, slice_type='x')
show_slices(x, y, z, volume_data, slice_type='y')
show_slices(x, y, z, volume_data, slice_type='z')

# convert to healpy format
update_antenna_conventions={
        'shift_phi':-90,
        'flip_theta':True,
        'flip_phi':False,
        'in_degrees':True,
        'add_invisible_sky':True
    }

antenna_hpmap_inst = antenna_inst.convert2hp(frequency=45, **update_antenna_conventions)

# hp maps are by default in galactic coordinates, so we need to specify the LST and latitude of the local observer
lst = 18
LATITUDE = -35.206667
rotation_parameters = create_rotation_parameters(lst, LATITUDE)
rotator = create_rotator(lst, LATITUDE, coord=["G", "C"])

antenna_hpmap = antenna_hpmap_inst.get_map(rotator=rotator)

# galaxy
projview(
    antenna_hpmap,
    cmap='jet',
    return_only_data=False,
    graticule=True,
    graticule_labels=True,
    title='Galactic coordinates',
    xtick_label_color='w',
    # projection_type='cart'
)

# from galaxy to local
projview(
    antenna_hpmap,
    cmap='jet',
    return_only_data=False,
    coord=['G','C'],
    rot=rotation_parameters,
    graticule=True,
    graticule_labels=True,
    title='Local coordinates',
    xtick_label_color='w',
    # projection_type='cart'
)